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PREFACE 


This  report  is  part  of  a  continuing  research  effort  sponsored  by 
the  Defense  Advanced  Research  Projects  Agency  (DARPA)  and  the  U.S. 

Arms  Control  and  Disarmament  Agency  (ACDA)  to  resolve  technical  issues 
concerning  verification  of  nuclear  test  ban  treaties.  Volume  I  of 
this  report  presents  a  procedure  for  estimating  station  seismicity, 
noise,  and  magnitude-bias  parameters.  The  noise  parameters  are  re“ 
quired  inputs  for  the  Seismic  Network  Assessment  Program  for  Detection 
(SNAP/D).  Volume  II  applies  the  procedure  developed  in  Vol.  I  to  the 
AEDS  classified  seismic  network. 

An  earlier  identically  titled  version  of  this  report  (PSR  Report 
1552,  Vols.  I  and  II,  August  1985)  derived  station  parameter  estimates 
without  accounting  for  the  effect  of  the  network  detection  criteria  on 
station  histograms.  This  version  accounts  for  this  effect  and  extends 
the  results  of  the  earlier  report  to  station-amplitude  histograms  as 
well  as  station  mjj  histograms. 
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I.  INTRODUCTION 


This  report  develops  a  model  that  can  be  used  in  conjunction  with 
single-station  histogram  data  to  estimate  station  seismicity  and 
performance  parameters.  The  histogram  data  generally  consists  of  a 
plot  of  the  number  of  earthquakes  (from  a  restricted  epicentral 
region)  versus  station  m^,,  although  a  model  for  the  treatment  of 
station  amplitude  histograms  for  worldwide  seismic  data  is  also 
considered.  The  estimated  station  parameters  consist  of  the  mean 
noise  level,  y,  the  standard  deviation  of  log  noise,  and  the 
station  magnitude  bias,  e.  In  order  to  estimate  e,  a  similar  his¬ 
togram  of  network  detection  performance  corrected  for  maximum- 
likelihood  (ML)  m^j  magnitudes  must  also  be  available. 

The  noise  parameters  y  and  0^  for  each  station  are  required 
inputs  for  the  Seismic  Network  Assessment  Program  for  Detection 
(SNAP/D)  [Ciervo,  et  al.,  19833.  Previously,  their  values  were  either 
estimated  from  measurements  made  on  seismograms,  or  inferred  from 
station  detection  thresholds.  The  empirical  method  does  not  generally 
ensure  accurate  replication  of  station  performance  in  SNAP/D  runs  and 
the  threshold-inferred  estimates  are  generally  reliable  only  for 
magnitudes  near  the  station  threshold.  The  procedure  presented  here 
is  based  on  past  station  performance  throughout  the  magnitude  range 
experienced  by  the  station  and  is  thus  believed  to  be  an  improvement 
on  prior  noise  parameter  estimation  procedures. 

A  similar  treatment  for  a  somewhat  different  problem  has  been 
presented  by  Kelly  and  Lacoss  [1969]  where  ML  estimates  were  derived 
for  network  performance  parameters.  However,  for  mathematical  con¬ 
venience,  a  single-station  Gaussian  detection  model  was  used  to  repre¬ 
sent  the  network  detection  process.  Furthermore,  the  biasing  effect 
of  non-ML  corrected  network  m^j  estimates  [see  Ringdal  1976]  was  not 
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understood  at  that  time.  A  similar  effect  on  single-station  seis¬ 
micity  estimates  is  accounted  for  in  the  procedure  derived  here.  In 
addition,  single-station  noise  estimates  are  corrected  for  the  usual 
four-station  P-wave  network  detection  criteria. 

The  procedures  developed  here  are  illustrated  using  histogram 
data  for  the  station  HYB  (Hyderabad,  India)  observing  events  from  the 
Kamchatka/Kurile  region  of  the  USSR.  The  estimation  procedures  are 
also  applied  to  the  analysis  of  classified  station  mjj  histograms  as 
detailed  in  Vol.  II  of  this  report.  However,  station  Mg  histogram 
data  was  either  too  sparse  or  irregular  to  obtain  reliable  estimates, 
hence  Vol.  II  presents  Mg  noise  parameters  from  a  previously  published 
report  [Hutchenson,  19833. 
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II.  HOTATION 

The  following  notation  will  be  adopted  for  the  discussion  below 

m  operational 

*■  observed  at  a  single  station 
a  -  log  amplitude  (log  A/T) 

a  »  log  amplitude  observed  at  a  single  station 

N(m)  -  density  for  the  expected  number  of  earthquakes 

occuring  with  magnitude  m 

a  »  intercept  of  base  e  seismicity 

B  ■  slope  of  base  e  seismicity 

p  -  station  mean  noise  amplitude 

e  «  station  magnitude  bias 

r  -  SNR  required  for  station  observation 

b(A)  »  b-f actor  (i.e.,  m  -  log  (A/T)  -  b(A)) 

p'  .  -  “b(A)  +  log  p  +-  log  r 

$  -  unit  normal  probability  distribution 

0  »  standard  deviation  (s.d.)  of  m  given  m 

s 

-  s.d.  of  single-station  log  noise 
yj^  ■  number  of  events  in  kth  magnitude  interval 
of  station  histogram 
i  -  station  index 

j  «  epicentral  region  index 
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III.  SEISMICITY  MODEL 


It  is  generally  accepted  that  logarithmic  seismicity  from  a  given 
region  Is  linear  with  respect  to  seismic  magnitude  [Richter,  1958], 
Thus,  if  N(m)A  is  the  average  number  of  seismic  events  per  year  occur- 
ing  in  the  operational  magnitude  interval  (m  -  A/2,  m  +  A/2),  then 
define  o  and  6  such  that 


N(m)  -  e“'"Bni  .  (1) 

As  in  Kelly  and  Lacoss  [1969],  the  actual  number  of  events  in  an 
operational  magnitude  interval  of  width  A  is  assumed  to  be  Poisson 
distributed  with  mean  N(m)A.  Although,  to  the  best  of  our  knowledge, 
no  formal  justification  has  been  offered  for  this  assumption,  it  is 
reasonable  since  over  a  fixed  time  interval  (say,  one  year),  the 
occurrence  of  primary  earthquakes  parameterized  on  magnitude  appear  to 
satisfy  the  axioms  of  a  nonhomogeneous  Poisson  process  [Parzen,  1962]. 

Define  the  kth  operational  magnitude  interval  as  (m)^  -  A/2,  m|^  + 

'  A/2),  vrtiere  A  «  mj^+i  -  mj^,  k  -  1,  2,  ...,  and  Xj^  as  the  random  number 
of  earthquakes  with  operational  magnitude  within  the  kth  interval. 

Then  the  Poisson  assumption  implies  that 

-N(m  )A  [N(m  )A]’' 

f{X.  -  x}  -  e  -  X  -  0,  1,  2,  ...  (2) 

K  XI 

where  N(mi^)  is  given  by  Eq.  (1). 

At  this  point  no  claim  has  been  made  about  the  distribution  of 
1)^,  the  number  of  earthquakes  detected  by  a  single  station  with  Ob'- 
served  magnitude  in  the  Interval  (%  “  ^^2,  +  A/2).  However,  using 

Eq.  (2)  and  the  results  below.  Appendix  A  proves  that  is  also  Poisson. 

*For  a  discussion  of  true,  operational,  and  observed 
magnitudes  see  von  Seggern  and  Blandford  [1976].  Unless  otherwise 
noted,  all  magnitudes  are  mij  values  with  the  subscript 
suppressed. 
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IV.  SINGLE^STATIOH  DETECTION 


The  amplitude  of  a  seismic  signal  arriving  at  a  station  may  be 
considered  to  result  from  a  series  of  random  multiplicative 
(attenuation)  effects  on  the  seismic  source  amplitude.  The  central 
limit  theorem  would  then  imply  that  the  log  of'  the  station  amplitude, 
and  hence  observed  magnitude  m,  is  a  Gaussian  randan  variable 
given  the  operational  magnitude  m.  Thus, 

,  /a  \2 

1  /  m-m\ 

1  "  ^\  a  / 

f  {m|m}  -  —1—  e  ®  '  (3) 

0  /2ir 

s 

where  O3  is  the  log  signal  s.d. 

The  log  of  the  station  noise  amplitude  at  any  time  is  also  Gaus¬ 
sian  since  the  noise  is  generally  composed  of  signals  from  myriad 
minor  seismic  disturbances.  Suppose  a  station  detects  a  seismic 
signal  with  amplitude  s  when  s/n  >  r  where  n  is  the  noise  amplitude 
and  r  is  the  signal-to-noise  ratio  (SNR)  required  for  detection.  The 
probability  of  detection  would  then  be 

Pd  -  ^{s/n  >  r} 

-  f{loz  s  -  (log  n  +  log  r)  >  0}  .  (if)* 

From  the  discussion  above,  log  s  and  log  n  are  Gaussian  with  expecta¬ 
tion  and  variance  given  by  (in  SNAP/D  notation  [Ciervo,  et  al.,  1983]) 

Edog  s)  -  log(A/T) 

V(log  s)  -  o| 

E(log  n)  -  log  y 


*The  following  notation  is  used:  log^o  •  log  and  logg  -  An. 
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and 


V(log  n)  =  . 


where  A  is  the  mean  signal  amplitude  in  nm  at  the  dominant  wave  period 
T.  Seismologists  prefer  to  use  the  queintity  A/T  because  of  its 
relationship  to  the  energy  in  the  wave  train  [Richter,  1958]. 

Thus,  Eq.  (4)  becomes 


$ 


log(A/T)  -  (log  u  -t-  log  r). 


2^2 
a  +  a 
s  n 


vrtiich  is  essentially  Eq.  (6)  in  the  SNAP/D  User's  Manual.  The 
relationship  between  magnitude  m  and  amplitude  A  is 


(5) 


m  -  log(A/T)  -  b(A)  (6) 

where  b  is  the  correction  factor  for  epicentral  distance  A.  Defining 


u’  »  -b(A)  +  log  y  +  log  r 


(7) 


Eqs.  (6)  and  (7)  allow  Eq.  (5)  to  be  rewritten  as 


P 


D 


(8) 


Note  that  Eq.  (8)  is  actually  the  probability  of  single-station 
detection  conditioned  on  operational  magnitude  m.  It  is  also  useful 
to  consider  the  detection  probability  conditioned  on  the  observed 
magnitude  m.  In  this  case,  the  only  uncertainty  is  the  noise 
amplitude  so  that  [Von  Seggerq  and  Blandford,  1976], 
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I  ffi}  -  $ 


where  ^  denotes  detection. 
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V.  STATION  HISTOGRAM  MODEL 


Incremental  histogram  data  is  generally  a  plot  of  y|^  -  the  number 
of  events  detected  with  observed  magnitude  in  the  interval  -  A/2, 

+  A/2)  \diere  A  -  k  -  1,  2,  ...  .  The  histogram  data  y;^ 

is  a  realization  of  the  random  variable  discussed  on  p.  4.  The 
expectation  of  Y)^,  which  is  needed  for  station  parameter  estimation, 
is  derived  below  from  the  seismicity  and  single-station  detection 
models  above. 

Recalling  that  N(m)  is  the  average  density  of  earthquakes  at 
operational  magnitude  m,  the  average  density  of  earthquakes  arriving 
at  a  station  with  observed  magnitude  fi  is 

^{m|m}  N(m)  dm 
o 


Thus,  the  average  density  of  earthquakes  detected  by  a  single  station 
is,  using  Eqs.  (1),  (3)»  and  (9),  given  by 


where  a*  «•  a  +  0.5  and  the  approximation  is  due  to  the  negligible 

effect  of  using  0  instead  of  -«  for  the  lower  limit  of  the  integral. 
Thus,  as  noted  by  von  Seggern  and  Blandford  [1976],  the  apparent 

-8- 


AC5VC1 1 0 


effect  of  using  station  histogram  data  to  estimate  seismicity  is  to 
introduce  an  upward  bias  of  0.5  into  the  intercept  parameter  a. 

The  expectation  of  is  given  by 


E(Y|^)  -  N(mk)  A 


(11) 


where,  from  Eq.  (10),  B  ••  ()i,n  10)”^  B  and  A'  ■  (Hn  10)  ^  (itnA  +  a'). 

The  unbiased  histogram  intercept  is  then  A  -  (in  10)“^  (in  A  +  a)  so 
that  the  upward  bias  is  A'  A  -  0.5  (in  10)“H2a|  -  0.5(in  10)B2a|. 
Using  the  assumption  that  operational  seismicity  is  Poisson,  Appendix 
A  proves  that  the  Y|^  are  Poisson  random  variables  so  that  Eq.  (11)  is 
also  the  variance  of  Yj^. 

The  above  discussion  assumes  that  detection  at  an  individual 
station  is  independent  of  the  performance  of  the  remaining  stations  in 
the  network.  In  fact,  the  teleseismic  detection  of  an  earthquake 
generally  requires  at  least  a  four'*station  detection  of  the  P~wave. 

The  effect  of  ignoring  network  influence  is  to  overestimate  station 
noise  levels.  To  properly  account  for  the  effect  of  the  network 
detection  criteria  on  the  station  i  histogram,  the  probability  of  the 
remaining  stations  in  the  network  (i.e.,  the  "reduced  network") 
providing  at  least  a  three-station  P-wave  detection  must  be  con¬ 
sidered.  With  this  approach,  Eq.  (10)  becomes 


N(m) 


|m}  N(m)  dm 


(12) 


where  ^3  and  21^  represent  detection  by  station  i  and  the  reduced 
network,  respectively.  Since  a  closed  form  expression  for 
does  not  exist,  SNAP/D  is  run  over  a  range  of  m's  for  each  reduced 
network.  The  integral  in  Eq.  12  is  then  computed  numerically. 

-9- 
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Note  that  each  SNAP/D  run  in  the  numerical  evaluation  of  Eq.  12 
requires  the  noise  inputs  u  and  0  for  each  station  in  the  reduced 
network.  This  requires  an  iterative  process  which  begins  with  zero^'^ 
order  p  and  0  estimates*  (from  the  assumption  of  independent  station 
histograms)  to  determine  ^ulm}  for  each  reduced  network,  which  in 
turn  are  used  to  calculate  first  order  y  and  0  estimates,  and  so  on. 
The  detailed  results  of  such  an  iterative  process  are  shown  in  Vol.  II 
of  this  report  for  the  AEDS  network.  After  five  iterations,  the  es¬ 
timated  values  for  y',  a,  and  A  converged  to  fifth  digit  accuracy. 

The  resulting  y  values  converged  to  at  least  three  digit  accuracy. 


The  parameter  estimation  procedure  is  described  in  Section  VII. 


-10- 
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VI.  AMPLITUDE  HISTOGRAMS 


The  expressions  derived  in  the  previous  section  can  be  general¬ 
ized  to  station  i  log  aimplitude  histograms  provided  the  event's 
epicentral  region,  J,  is  known,  along  with  the  seismicity  of  the 
region.  Denoting  the  quantity,  log  A/T,  as  'a'  for  brevity,  Eq.  (6) 
yields 

a  »  m  +  b(A) 


The  attenuation  factor,  b,  depends  on  the  angular  distance  between 
station  i  and  epicenter  j,  which  is  indicated  in  the  following 
development  by  b^j .  Noting  that  the  observed  m^  is  related  to  the 
operational  m^  by  a  normally  distributed  deviation  with  mean  e^j 

and  standard  deviation  03^,  the  above  expression  becomes 


*i  ■  "1  *  •’ij "  • 


Eq.  (12)  can  be  now  be  rewritten  as 


where 


r  N.(m)  Jj,|m  at  j}  |m  at  J)  (13) 

j  L  ■ 


o  +B  m 

Nj(m)  -  e  ^ 


-1/2 

^{a.  |m  at  j}  -  ® 

/2ir  0  , 
si 


\  ~  ”11  ~  °ii  ~  ^ 
°si 


and 


C-^) 
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The  necessary  inputs  for  such  a  general  expression  would  be  a j ,  $ j , 
ej[j,  and  Ogj.  The  integral  still  must  be  evaluated  numerically 
with  computed  by  SNAP/D. 

There  are  three  special  cases  for  the  amplitude  histogram 
problem: 

1.  »  B  for  all  j:  with  this  aissumption,  B  can  be  estimated  as  in 
the  earlier  discussion  by  treating  data  as  independent  and 
averaging  the  fitted  Bj's.  The  SNAP/D  iterations  and  the 
numerical  integral  would  still  be  required. 

2.  Independent  amplitude  histograms  are  equivalent  to 

jj|ni  at  j }  H  1.0 

and  would  remove  the  need  for  the  SNAP/D  iterations. 

3.  Independent  amplitude  histograms  with  Bj*  B  would  cause  Eq.  13 
to  become 


Vfhere 

fl 

-  log  +  log  r^^ 


Defining 


“j  ^ 


2 

®  "si 


B  e 


ij 


1  2- 
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so  that  the  functional  form  of  the  expected  number  of  events  in  a  bin 
is  similar  to  the  independent  m^  case  Eq.  (10).  In  Eq.  (10),  o',  yj, 
and  aj  would  generally  be  regarded  as  free  parameters  to  be  estimated, 
while  6  may  be  regarded  as  fixed  or  free. 
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VII.  PARAMETER  ESTIMATION 


For  ease  of  computation,  a  minimum  chi-square  (MCS)  estimation 
procedure  was  chosen.  However,  Appendix  B  proves  the  asymptotic 
(large  I-y|<)  equivalence  of  MCS  and  ML  estimates  for  this  problem. 
For  each  station  the  MCS  estimates  are  those  that  minimize  the  usual 
chi-square  sum,  i.e.. 


(16) 


where,  in  our  case,  the  kth  observation  oj^  -  yi<,  and  the  expected 
observation  e|^  ■  N(m|^),  where  N  is  given  by  Eq.  (10)  or  (12)  depending 
on  whether  independent  or  dependent  station  data  is  used  .  In  prac¬ 
tice,  the  MCS  sum  is  computed  over  all  k  for  which  ^  1 ,  where  the 
parameter  values  at  each  stage  of  the  iterative  minimization  process 
are  used  to  compute  ej^.  The  technical  properties  of  the  MCS  estimates 
(and  their  asymptotic  equivalence  to  maximum  likelihood  (ML)  es¬ 
timates)  are  detailed  in  Appendix  B  where  the  variance  of  the  es¬ 
timates  are  also  derived. 

Estimates  for  each  station  are  calculated  in  three  steps:  (1)  the 
minimization  indicated  in  Eq.  (16)  is  performed  (under  the  assumption 
-  of  Independent  station  histograms)  with  all  four  parameters  free;  (2) 
the  weighted  average 


i 

is  calculated  where  Bj^  is  the  station  i  slope  estimate  from  step  (1), 
nj  is  the  number  of  events  iq  the  station  i  histogram  and  n  -  Zn^;  and 
(3)  the  minimization  is  repeated  with  three  free  parameters  and  B  »  F 
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fixed  for  all  stations.  In  the  case  of  network  dependent  histograms, 
the  iterative  minimization  procedure  presented  in  Sec.  V  must  be 
used.  The  rationale  for  setting  B  »  ?  to  obtain  the  final  u' ,  and 
A'  estimates  is  that  B  is  generally  a  better  estitnate  of  the  opera¬ 
tional  slope  than  an  individual  Bj. 

SNAP/D  requires  mean  station  noise  levels  in  amplitude  units  (0-P 
in  nm)  and  standard  deviations  in  log  amplitude  units.  Thus,  if  the 
minimization  procedure  described  above  provides  and  estimates 
for  station  i,  then  the  SNAP/D  station  i  noise  inputs  are,  from 
Eq.  (7), 


Mi  -  10 

and  Ojji  unchanged. 

The  complications  arising  from  network  influences  on  m^  his¬ 
tograms,  as  discussed  in  Sec.  V,  imply  that,  in  this  case,  an  unbiased 
seismicity  estimate  may  not  be  available  from  the  MCS  fit.  However, 
the  results  of  applying  the  procedure  to  AEDS  station  data  (presented 
in  Vol.  II  of  this  report)  indicate  a  negligible  change  in  A* 
parameter  estimates.  Thus,  the  expression  for  the  unbiased  station  i 
intercept  parameter,  Aj  ■«  A’j  **  O.SCAn  10)B^o^,  as  discussed  in  Sec.  V 
for  the  case  of  independent  m^  data,  is  also  used  as  a  reasonable 
approximation  for  the  case  of  dependent  data.  Since  SNAP/D  calcula¬ 
tions  are  conditioned  on  the  occurrence  of  a  seismic  event,  seismicity 
is  not  a  SNAP/D  input  for  unassociated  data,  but  estimation  of  the 
intercept  parameter  permits  estimation  of  station  magnitude  bias  as 
discussed  below. 

If  network  estimates  of  operational  seismicity  parameters  Ajjgj 
and  Bjjet  are  available,  then  the  MCS  procedure  would  consist  only  of 
step  (3)  with  Bj^  set  equal  to  Bjigf.  In  this  case,  the  estimate  of 
station  magnitude  bias  would  be 

’  %ET  -  *i 

“  Bjjet 
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this  expression  assumes  that  the  station  i  was  in  operation 
continuously.  If  the  periods  for  which  station  i  is  inoperable  are 
known,  an  obvious  modification  of  Ajjg-p  could  be  made  so  that  the  e^ 
calculation  would  still  be  correct.  The  usual  caution  associated  with 
using  network  histogram  data  to  estimate  seismicity  must  be  observed: 
only  MLE  corrected  network  magnitude  data  can  be  used  to  obtain  un- 
biased  estimates  of  Ajjgj  and  due  to  the  "bulge"  phenomena  as¬ 

sociated  with  histograms  based  on  network  average  magnitudes  [Zavadil, 
et  al.,  1983].  Although  SNAP/D  can  in  principle  accommodate  correc¬ 
tions  for  station  magnitude  bias  through  use  of  the  correction  factor 
epji^  (in  SNAP/D  notation:  i  -  station  index,  j  =  epicenter  index,  and 
k  •  wave  index),  in  practice  the  required  data  acquisition  and 
analysis  would  be  formidable  even  when  restricted  to  P-wave  observa¬ 
tion  for  seismically  active  areas  in  the  Soviet  Union. 
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VIII.  EXAMPLE 


The  calibration  procedure  described  above  was  applied  to  1976- 
1980  histogram  data  compiled  by  F.  Ringdal  [Rivers,  1984]  for  Kam¬ 
chatka  events  observed  by  the  station  in  Hyderabad,  India.  The  ap¬ 
plication  assumes  independent  station  data  for  simplicity  of  discus¬ 
sion.  This  is  clearly  an  approximation  to  the  true  case  of  network 
association.  Figure  1  plots  the  data,  the  MCS  fit  (Eq.  (6)),  the 
unbiased  station  seismicity  ^nd,  since  no  ML  network  mag¬ 

nitude  data  was  available,  hypothetical  unbiased  network  seismicity 
(|oAnet'^%ET'").  Assuming  B^et  -1,  A  -  0.1,  and  O3  -  0.35,  the  sta¬ 
tion  magnitude  bias  0.5(iln  10)b2o|  is  also  indicated.  Note  that  the 
apparent  station  seismicity  bias  refers  to  the  "vertical"  difference 
between  apparent  asymptotic  (large  magnitude)  station  seismicity 
and  the  unbiased  station  seismicity  On  the 

other  hand,  the  station  magnitude  bias  refers  to  the  "horizontal" 
difference  between  the  unbiased  station  seismicity  and  the  unbiased 
network  seismicity  (io^NET'*‘%ET®) . 

The  MCS  estimates  for  HYB  are 

M'  -  4.83 
“  0 .20 
A*  -  6.15 


with  seismicity  slope  fixed  at  B  -  -0.89.  Thus,  SNAP/D  noise 
parameters  would  be 

y  -  10l^*  +  b(A)  -  log  r 
-  7. ,48 
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Figure  1.  Illustration  of  calibration  fit  and  parameters. 
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and  On  -  0.20  where  b(A)  -  153.78  is  the  Veith  and  Clawson  [1972]  0  km 
entry  for  A  -  68.82°  and  r  -  1.5.  The  unbiased  seismicity  intercept 
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Appendix  A 

THE  DISTRIBUTION  OF  OBSERVED  MAGNITUDES 


This  appendix  will  show  that  if  X(m),  the  random  number  of 
earthquakes  with  operational  magnitude  in  the  interval  [0,m),  is  a 
Poisson  process,  then  Yjj(ffi),  the  number  of  earthquakes  detected  by  a 
single  station  with  observed  magnitude  in  [0,(9),  is  also  Poisson. 
Hence,  the  histogram  data  Y|^  discussed  on  p.  (t  of  the  text  is  a  Pois¬ 
son  random  variable. 

Let  X(m)  be  a  Poisson  process  with  density  function 

\  f  _ct— 6m 
Atm;  “  e  , 


where  6  >  0,me[0,*),  and  6  «  -S,  in  the  notation  of  Sec.  III.  X(m) 
denotes  the  total  number  of  earthquakes  with  operational  magnitude  in 
the  interval  C0,m),  and  X(m)  is  the  total  number  of  earthquakes  with 
magnitude  bigger  than  m  (i.e.,  magnitude  in  the  interval  (m,*))  .  Both 
X(m)  and  X(m)  are  assumed  to  be  Poisson  processes,  with  respective 
mean  value  functions 


E[X(m)]  -  A(m) 


m 


A 

^ _ 

X(m)  is  the  Poisson  process  treated  by  Kelly  and  Lacoss  [1969]. 
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ECX(m)]  -  A(m) 


0» 

/ 


A(s)  ds 


,o-6m 


Is 


In  an  interval  ^ni  -  m  +  the  expected  number  of  earthquakes 


-  -a“6ra  . 

*  e  A 

Now  to  develop  the  distribution  of  Y(ffi),  the  number  of 
earthquakes  with  observed  magnitude  in  [O.fil).  let  be  the  opera¬ 
tional  magnitude  of  the  i^'^  smallest  earthquake  recorded  by  the  sta¬ 
tion,  where  may  be  smaller  or  bigger  than  fl.  Let  denote  the 

corresponding  observed  magnitude,  so  that 


fii  ■  hi  +  Ei  , 


where  the  measurement  error  is  Gaussian  with  mean  0  and  standard 
deviation  Og. 

Define  the  0-1  valued  function 


w(l!l,n,e) 


(  1  if  n  +  e  S  la 
\o  if  n  +  £>fl 


Then  the  i  earthquake  is  included  in  the  Y(fi)  count  if  fij  i  ifi,  i.e., 
if  w(fli,  T|j,  »  1,  and  otherwise  not.  Thus,  Y(fil)  may  be  expressed 
as  « 
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Y(ffi)  - 


Clearly,  Y(ft)  is  a  nonnegative  integer-valued  random  variable. 

The  assertion  that  Y(ft)  is  again  a  Poisson  process  can  be  proved 
in  several  different  ways,  but  we  shall  present  a  proof  that  is 
straightforward  and  intuitively  appealing. 

First,  we  demonstrate  a  fairly  well  known  lemma  [Papoulis,  1984] 
relating  the  limiting  distribution  of  a  sum  of  independent  Bernoulli 
(0-1  valued)  random  variables  to  a  Poisson  variate. 

Let  x^,  X2,  ...,  x^,  Xj^+i . be  a  sequence  of  independent 

random  variables  such  that  x^^  is  1  with  probability  p^  and  0  with 
probability  -  1  -  p^.  Further,  let 

n 

X  -  lim  7^  Pi  <  * 

n^co 

and  assume  that  max  Pi->-0  as  n-*-". 

12l^n 


Then  we  assert  that  Zj^ 


converges 


(as  iiH-oo)  to  a  Poisson  ran¬ 


dom  variable  Z  with  mean  X. 

To  show  this,  we  exhibit  the  characteristic  function  of  each 
and  of  Z,  and  show  consequently  that  the  characteristic  function  of  Z„ 
converges  to  that  of  Z,  thus  establishing  the  lemma. 
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The  characteristic  function  of  the  Bernoulli  variable  Xj  is 


♦,(u)  . 


(where  J 
Z  is 


-  Pi  eJ’^  +  Qi 

and  the  characteristic  function  of  the  Poisson  variate 


(|)2(u)  - 


Juk 


k-0 


k! 


k-0 


k! 


,A(eJ“-1) 


Now  note  that  if  <<  1 ,  then 


j  ^  s  1  +  Pi(eJ’^-l) 


PieJ“  +  Qi 


-  <|>i(u) 


(A.1) 


and  so 
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n 

£n>2„(u)  -  V] 

In  (frj^(u) 

(by  independence  of  Xj) 

i-1 

n 

■E 

4n(pj^e'^“  +  1 

■  Pi) 

i-1 

Hr-  -1 

■E 

p.p-') 

1  +  e^Pi 

i-1 

L 

J 

n 

n 

■  E 

Pi  (eJ“  -  , 

i)  *  2  *1  Pi 

i-1 

i-1 

-  1) 

(A. 2) 

where  6j;+0  as  Pj->-0.  Hence  (j>2  (u)— ►<{>2(u) ,  establishing  the  lemma. 

To  show  that  if  X(m)  is  Poisson,  then  so  is  Y(fi),  begin  by  par¬ 
titioning  the  magnitude  axis  C0,“)  into  consecutive  intervals  ■ 

Oi+i)  of  fixed  length  (or  mesh)  Ao  -  -  ot^,  as  in  Fig.  A.1. 

Let  AXj^  denote  the  number  of  earthquakes  with  operational  magnitude  in 
the  interval  and  let  AY (ft,  oj)  denote  the  corresponding  contribu¬ 
tion  to  the  sum  Y(ft)  due  to  the  error  in  measuring  the  operational 
magnitudes  of  the  earthquakes  in  interval  1^.  Thus,  we  may  write  Y(ft) 
as 


Y(ft) 


AY(ft,  o^) 
i 


0 

«2 


i — I — 


Figure  A.1 
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If  the  mesh  Ao  of  the  partition  is  small,  then,  within  probabil¬ 
ities  of  order  La,  the  random  variable  AXj  takes  the  value  0  or  1 ,  and 

P(AXj  -  1)  =  XCoj)  Ao 

where,  as  before,  X(m)  -  e®”*®  is  the  intensity  function  of  the  Pois¬ 
son  process  X(m) . 

If  LX^  -  0,  then  necessarily  AY(fil,  a^)  *  0.  But  if  *  1 ,  then 
AY(ft,  a^)  *  1  <r — >  ^  fi 

or  ^  fi  -  Tlj^^ 

where  the  operational  magnitude  lies  ii^  interval  Ij^. 

Thus,  approximately,  given  that  AXj^  •  1 , 


p[AY(ia,ai)  -  1  1  AXi  -  lj  =«!• 

since  the  measurement  error  is  Gaussian  with  mean  0  and  standard 
deviation  o_.  Hence  the  unconditional  probability  is 

O 


•[^AY(fl,  Oi)  “  l]  =  ^(“i^ 


Now,  the  random  variables  AY(fi,  aj^)  take  the  values  0  or  1  and 
are  independent,  because  the  are  independent  (Poisson  variates  in 
nonoverlapping  intervals)  and  the  measurement  errors  ^  are 
independent.  Therefore,  by  the  lemma,  as  the  mesh  Aa-^,  the  sum  Y(fi) 
-  I  AY  (ft,  a^)  tends  to  a  Poisson  variate  with  mean 


lim  XCa^)  ^(7 - ^  La 

Ao-»0  \  ®s  / 


(A. 3) 
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The  above  argument  also  illustrates  that  for  two  nonoverlapping 
intervals  (0,  S)  and  (fi,  A  +  t),  where  S  <  A,  the  variates  Y(S)  and 
Y(A  +  t)  -  Y(A)  are  independent,  and,  of  course,  Poisson.  Hence  Y(A) 
is  also  a  Poisson  process,  as  was  to  be  shown. 

From  Eq.  (A. 3),  we  may  write  the  mean  value  function  for  Y(A)  as 


and  the  corresponding  intensity  function  is 


(A-m)^ 

u(A)  »  —  M(A)  »  f  — =  e  ®  dm 

dA  J 

o  ® 

(A.4) 

N(m)  ?‘{A|m}  dm  , 
as  in  Section  V,  pg.  7. 

Finally,  we  consider  the  effect  of  the  probability  of  detection 

(A  ~  u' 

— 

on  the  process  Y(A).  As  before,  partition  the  magnitude  axis  into 
disjoint  intervals  ■  (Sj^,  S^+i )  with  mesh  AS  ■  5^+^  -  Sj^,  and  let 
AS  be  so  small  that  within  probabilities  of  order  AS  the  number  AY^  of 
occurrences  of  earthquakes  with  measured  magnitude  in  is  0  or  1 . 
If  AYj  -  0,  define  AY^CA,  S^)  -  0. 

If  AYj  -  1,  define 

AYjj(A,  Sj)  -1*  < — ►  the  event  in  1^  is  detected. 
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Thus,  given  AYj  -  1,  the  conditional  probability  that  AY^Cfil,  is  1 
is 


pjAYu  (fi,  a^)  -  1  I  AYj  "  ■'j  ‘ 

and  as  developed  earlier 


P(AYi  *  1)  ■  p(a<)  Ao  . 


The  total  number  of  earthquakes  detected  in  [0,  ft)  is 


YoCft)  -  ^  AYijCft,  a^)  . 


And  from  Eqs.  (A. 5)  and  (A. 6), 


'jAYjj(ft,  a^)  *  “  u(aj^)^{^|aj^}  Aa^ 


(A. 5) 


(A.6) 


Thus,  as  the  mesh  Aa^0,the  total  number  of  earthquakes  with  magnitude 
in  [0,  ft)  that  are  detected  is  a  Poisson  process  with  mean 


(A. 7) 


Mjj(ft)  -  f  u(x)  dx 

•b 

In  an  interval  ^ft  “  ®  number  of  earthquakes  that  are 


detected  is  given  by 


%(*  *  I)  -  »d(*  -  1)  -  /.  y(x)  dx 

ft-^ 
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f  y(x)  dx 

4-^ 

Finally,  by  differentiating  Eq.  (A. 7),  we  obtain  the  density  (see 
Eq.  (A. 4)  for  iJi(fl))  for  the  detected  process, 


UpCfa)  -  y(fi)  fVjD\^) 


N(m)  dm 


“  N(fil)  , 

as  in  Eq.  (10)  of  Sec.  V  of  the  text. 
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Appendix  B 

MINIMUM  CHI-SQUARE  AND  MAXIMUM  LIKELIHOOD 
IN  FITTING  A  POISSON  PROCESS  MODEL 


In  this  appendix,  we  assume  that  count  observations  can  fall  into 
any  one  of  K  bins  or  magnitude  intervals,  where  each  bin  has  the  same 
width  A.  Let 

yj^  =  number  of  observed  counts  in  bin 
ej^  «  expected  number  of  counts  in  k*^  bin 

/"k  -  “X 

=  e“  h(m|^,  0') 

=  Yh(m^,  0) 

for  each  k,  k  =  1,  ...,  K,  where  m^^  is  the  midpoint  of  the  kth  bin, 
0  =  (0,  u,  Cjj),  0’  *  (o,  0),  and  the  definition  of  h(.)  is  clear  from 
the  above.  Note  Y  =  e“. 

Assuming  that  y^  yj^  are  independent  Poisson  random  vari¬ 

ables,  where  y|^  has  mean  A|^,  then  N  =  lyi^  is  also  Poisson  with  mean 
Moreover,  the  joint  density  function  (or  likelihood  function)  of 
y^,  ...,  yjf  is 


L(0’) 


f*(yi,  ....  yK) 


n  -X 
k  e 


k,  yk 


n  yu! 
k  ^ 


The  density  function  for  the  sum  N  is 


(B.1) 
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P(N  -  n) 


for  n  =  0 ,  1 ,  ...  . 


The  log  of  the  likelihood  function  from  Eq.  (B.1)  is  then 


log  L  -  -Z  Xk  *  Z  an  A;,  -  Z  4n  y^! 


Z  Yh(mj^  0)  +  Z  y,^  [an  Y  +  an  h(m,^,  0)1  -  Z  an  y^! 

k  k  t  J  k 

(B.2) 


MAXIMUM  LIKELIHOOD  ESTIMATION  (MLE) 

The  maximum  likelihood  estimate  (MLE)  for  0'  is  obtained  by 
choosing  to  maximize  L(0’)t  or  equivalently  log  L(0’).  If  L(0')  is 
differentiable  with  respect  to  0’  (as  in  our  model),  the  MLE  §'  is 
obtained  by  solving  the  equation 


d 

-  log  L(0’)  =•  0  .  (B.3) 

30’ 


By  differentiating  Eq.  (B.2)  with  respect  to  Y,  we  obtain 


or 


0 


3  log  L 
3Y 


1 

-  Z  h(m^,  0)  +  “  Z  yi.  , 
k  Y  k  '' 


Y  -  -  - - 

Z  h(m|^,  0)  Z  h(m|^,  0) 


(B.4) 


Substituting  this  expression  for  Y  in  the  log  likelihood  (B.2),  we 
obtain  (aside  from  a  constant  depending  on  N,  but  not  depending  on  0) 
the  log  likelihood  for  the  conditional  distribution  of  (y-j,  ...  ,  yj^) 
given  N  -  Z  y|^,  which  is  multinomial  with  parameters  N,  ir-|(0),  ..., 
■wk(0),  where  for  1  S  k  ^  K, 
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h(m,,,  0)  Aw(e') 

irv(0)  -  - - -  -  - 

E  h(m.,  6)  E  X.(e') 

J  J  J  J 

This  version  of  the  log  likelihood  may  then  be  maximized  (by  solving 
for  3  log  L/30  =  0)  to  obtain  the  MLE  estimates  for  0  =  (B,  u,  a^j). 

MINIMUM  CHI-SQUARE  ESTIMATION  (MCSE) 

Although  the  asymptotic  properties  of  MCS  are  well  known  for  the 
multinomial  case  [Cox  and  Hinkley,  1977  and  Rao,  1957],  they  do  not 
seem  to  be  not  as  well  known  for  the  Poisson  case.  Since  we  were  not 
able  to  find  in  the  literature  an  exact  reference  for  the  Poisson 
case,  we  present  the  development  here. 

For  the  Poisson  model,  the  MCS  estimate  is  achieved  by  minimizing 


S(0')  -  ^  ■“ 

^  ®k 

(yk  ~  X|^(9*))^ 

^  (B.5) 

where  ej^  is  the  expected  number  of  counts  in  bin  k,  which  in  this 
model  is  X|^(0’). 

The  MCSE  is  obtained  by  solving  the  equation 


3 

-  S(0’)  -  0  . 

30' 


We  will  show  that  under  suitable  regularity  conditions,  the  solution 
of 


3 

—  S(Y)  -  0 
3Y 


yields  an  estimate  Y  which  is  asymptotically  equivalent  to  the  ML 
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estimate  Y,  and  that  substitution  of  the  resulting  estimate  into 
Eq.  (B.5)  yields  (approximately)  the  MCSE  for  the  multinomial  dis¬ 
tribution,  which  is  known  to  be  asymptotically  equivalent  to  the  MLE. 
Now  the  MCS  criteria  is 


S(e') 


Lyk  ^ 

Yh(mi^,  e) 


(B.6) 


Letting  h|^  -  h(m|^,  0)  and  differentiating  with  respect  to  Y,  we  obtain 


hence 


is  the  MCS  estimate  of  Y,  for  fixed  0. 

The  ML  estimate  9  is  unbiased,  since 


Its  variance  (since  N  is  Poisson)  is 


t 
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The  variance  becomes  small  if  E  is  large,  which  we  shall  assume. 

k 

To  evaluate  the  mean  of  the  MSCE  y,  let 


y  Ik 

I,  h,, 

Z  -  — - ^ 


E  h. 

j  ^ 


Since  is  Poisson  with  parameter  Y  h,^, 

E  (y^)  -  Y^  h^  +  Y  h,^  . 

Thus 


I 


E(Z)  -  ^ 

E  h. 

J  J 


To  first  order,  then,  since  Y  = 
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where  the  last  approximation  follows  from  the  Taylor  expansion  of 


R.  _1  . 

If  -  goes  to  0  (or  K  Z  hi.  grows  large)  the  bias  term  goes 

Z  h,,  k 

k 


to  zero. 

It  may  also  be  shown,  using  the  first  four  moments  of  the  Poisson 
distribution,  that  to  first  order. 


Var  Y  s 


Z  h. 


just  as  for  the  MLE  Y. 

Thus,  if  K  and  Z  h.  grow  large  in  such  a  way  that  y/t  h, 
k  ^  k  ^ 

and  h^^  decreases  toward  zero,  V  and  Y  will  be  asymptotically 
equivalent. 


Writing 


Y  -  Y  +  e(Y) 


N 

Z  h. 


+  e(Y) 


where  the  error  e(Y)  depends  on  Y  and  on  Z  h|^,  we  replace  Y  by 
Y  -  N/Z  h^  +  eCY)  in  Eq.  (B.6),  and  obtain 
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But  if  |e(”lf)  h|^|  is  small  for  each  k,  this  yields  (approximately)  the 
MCS  criterion  for  the  multinomial  distribution  of  (yi  yj^)  given 
N.  From  Ferguson  [1958],  the  resulting  MCS  estimate  for  0  «  (B,  u,  a^j) 
is  asymptotically  equivalent  (as  N  grows  large)  to  the  ML  estimate 
(i.e.,  8  is  Best  Asymptotically  Normal,  or  BAN). 

The  asymptotic  equivalence  of  MCS  and  ML  estimates  allow  us  to 
use  the  usual  MSE  approach  [Cox  and  Hinkley,  197^]  to  evaluate  the 
asymptotic  variance  for  the  MCS  estimates.  Thus,  we  evaluated 


log  L(e’) 

0'  -  0' 

where  8*  is  the  MCSE,  and  inverted  this  matrix  to  achieve  the  es¬ 
timated  variance-covariance  matrix  for  0'.  Elements  of  this  matrix  are 
given  below  (with  a  sign  reversal): 


3^(logL(8* )) 
3a  3y 


♦ 
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In  summary,  we  have  shown  that 

1.  The  MLE  Y  is  unbiased  for  Y,  with  variance  inversely  propor¬ 
tional  to  Z  h|^;  thus  Y  converges  to  Y  in  probability  if  Z  hj^ 
grows  large. 

2.  The  MCSE  Y  is  biased,  with  bias 

K 


b(e) 


2  Z  h, 


As  K  ^  Z  hi.  grows  large,  the  bias  decreases  toward  zero.  The 
variance  of  Y  is,  to  first  order  Y/Z  hj^,  just  as  for  the 
MLE. 

We  may  write  y  -  Y  +  e(Y),  where  the  error  term  e(Y)  is  the 
sum  of  two  components 


i)  the  bias  term  b(0) 


K 


2  Z  h. 


ii)  a  mean  zero  random  component  with  variance 


Z  h. 


Asymptotically  as  Z  hj^  grows  large  and  Z  \/K  grows  large, 
the  error  e(Y)  becomes  small. 

Replacing  Y  by  Y  in  the  equation  for  the  MCSE  criterion 
yields,  because  of  continuity,  a  solution  that  is  asymptoti¬ 
cally  equivalent  to  the  MLE  for  6  =•  (B,  y,  a^j). 
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